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We develop an algorithm that has the potential to relate the depth development of ultra high en- 
ergy extensive air showers and the time delay for individual muons. The time distributions sampled 
at different positions at ground level by a large air shower array are converted into distributions 
of production distances using an approximate relation between production distance, transverse dis- 
tance and time delay. The method is naturally restricted to inclined showers where muons dominate 
the signal at ground level but could be extended to vertical showers provided that the detectors used 
can separate the muon signal from electrons and photons. We explore the accuracy and practical 
uncertainties involved in the proposed method. For practical purposes only the muons that fall 
outside the central region of the shower can be used, and we establish cuts in transverse distance. 
^■"^ . The method is tested using simulated showers by comparing the production distance distributions 

' obtained using the method with the actual distances in the simulated showers. It could be applied 

^ in the search for neutrinos to increase the acceptance to highly penetrating particles, as well as 

for unraveling the relative compositions of protons and heavy nuclei. We also illustrate that the 
' obtained depth distributions have minimum width when both the arrival direction and the core 

position are well reconstructed. 
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' When an Ultra High Energy Cosmic Ray (UHECR) particle enters the atmosphere it interacts producing an 
, extensive air shower that propagates through it and reaches ground level. These showers are routinely detected by 
' optical systems that collect fluorescence light emitted by nitrogen molecules excited as the front crosses the atmosphere, 
, and by arrays of particle detectors that sample at ground level an enormous shower front which can exceed 10^^ 
■ particles. In these arrays the relative times of the detected signals allow the reconstruction of the incoming particle 
[ arrival direction. The time distribution of the arriving signal has been known for long to be dependent on the depth 
O • distribution of the shower particles 1, which is different for different primary particles. 

I ' Exploring the highest energy particles is now considered to be a priority because the data are scarce, there are 
discrepancies between results obtained with the two techniques and because the origin of these particles is not at 
all understood (5|. Their study is expected to provide both information on violent objects in the Universe where 
^ , these particles originate and on their interactions (during propagation and in the atmosphere) at energies exceeding 
those achieved in accelerator experiments by many orders of magnitude. New experiments are being constructed and 
devised to improve the statistics, to increase the precision and to establish the mass of the primary particles. The 
Auger Observatory in Argentina is the first of a new generation of large aperture experiments. It combines the two 
5— i techniques and for the ground array it uses water Cerenkov tanks with photomultiplier tubes and Flash Analog to 
Digital Converters (FADC) to record the time stamp of the signal in 25 ns intervals with unprecedented accuracy Q. 

The perspective of improved detectors has triggered an increase of phenomenological activity in the study and 
characterization of extensive air showers. Evaluation of the time structure of showers is part of this effort motivated 
by both the practical need of controlling the uncertainties in the arrival direction reconstruction and also by the 
hope that its understanding might shed new light on the challenging problem of establishing their composition. The 
idea of relating the muon distributions to the shower development has been already quite successful 0, H, S ^3 ■ 
The arrival time distributions of muons has been characterized using simple geometrical and kinematical arguments 
and making a key simplifying assumption on the muon energy, transverse momentum, and distance of production 
distributions, namely that they are independent 12,]. Here we have further developed the algorithm that relates the 
arrival time distribution of muons to the depth development of the shower following on these ideas. 

The scope of the method is limited because the shower front contains many other particles, mainly electrons and 
photons. In principle it requires muon identification but fortunately the muon signal dominates in two circumstances: 
when the shower is inclined and the electromagnetic part does not reach ground level ^| but also for "more vertical" 
showers when the distance to shower axis is sufficiently large [T^ . The method complements alternative depth 
reconstruction methods which are always limited when only densities at ground level are taken into consideration. 
Moreover when the effects of arrival angle and impact point misreconstruction are taken into account it is seen 
that the induced distributions have a minimum in their spread for the correct angle and impact point. This effect 
opens the possibility of using this method for improving the confidence in the conventional angular and impact point 
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reconstructions. While the precision obtained is possibly insufficient to be used for composition studies it will certainly 
have an important impact on improving the acceptance of air shower arrays to neutrinos through inclined showers. 
The accuracy in the depth development reconstruction is sufficient to exclude neutrino interactions at intermediate 
depths, when the electromagnetic shower would have been completely absorbed but the first interaction is sufficiently 
deep into the atmosphere to exclude both a cosmic ray hadron and a photon. 

The article is organized as follows: In Section II we summarize the factorization hypothesis for the muon distributions 
and the relations that follow, and motivate the inversion of the relation between the time and depth distribution from 
Ref. [12]. We also pay some attention to the relation between particle densities and detector geometric acceptance. 
In Section III we present the method to reconstruct the depth distribution. In Section IV we compare the depth 
distributions obtained with this reconstruction method to actual distributions from simulated showers to test it. In 
Section V discussing some practical limitations. In Section VI we explore the correlations between the reconstruction 
procedure and the assumed arrival direction and impact point as a check of its robustness; we summarize and conclude 
in Section VII. Technical details are presented in two appendices. 



II. RELATION BETWEEN DEPTH DEVELOPMENT AND MUON TIME DISTRIBUTIONS 

The main features of the arrival time distributions of muons in extensive air showers can be accounted for by 
the different path lengths traveled by the muons from their production point. This has been recently shown using a 
simplified model to describe the muons in air showers which is based on the hypothesis that their energy, Ei , transverse 
momentum, pt, production distance, z, and outgoing polar angle, C distributions factorize |l2j| 

=^^foh{z)h{E,)f2{pt), (1) 



dz dC dEi dpt 271 

In this expression Ei^ pt and Aq refer to production, z is measured along the shower axis from the muon production 
point to the ground. The transverse momentum is transverse to shower axis and has polar angle C in the transverse 
plane (perpendicular to shower axis). The functions h{z), fi{Ei) and f2{pt) are assumed independent and normalized 
to 1 and the factor 27r accounts for a uniform polar angle distribution. Finally Ao, the total number of produced 
muons, is the overall normalization. In this model the muons are assumed to travel in straight lines and these four 
variables are sufficient to determine the muon path uniquely. 

It is convenient to express the muon direction in terms of the angle a with respect to shower axis. For a muon 
produced with energy E^ and transverse momentum pt, the angle a is given by 

Ptc _,PtC ^2) 
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We can approximate pc ~ ^ Ef — (mc^)^ ~ Ei, because the muon energy at ground level is typically greater than 
mc^, this energy at production is even greater because of muon energy loss. We can now change the coordinates 
replacing pt in Eq. ^ using Eq. [21 to give 

=Mo^hiz)ME,)f2(^sina)^. (3) 



dz dC dEi dsin a 2n \ c J c 

As the muons go through the atmosphere they lose energy and decay and even though we start with independent dis- 
tributions at production, correlations between the relevant variables appear naturally when we consider the surviving 
muon distribution at ground level. This is explicitly shown in Appendix A using a simplified model for energy loss. 

The distribution of surviving muons can then be integrated in Ei in order to obtain the depth distribution of the 
surviving muons, which is given in terms of two angles describing the muon direction, namely a and C. It is convenient 
to relate them to the differential solid angle for the muon d^il = —d(d cos a. Then 

d^N d^N cos a 

(4) 



dzd'^fl dzd(^dsina sina 

It is interesting to discuss in some detail the effect of a detector surface. From Eq. 2]we can obtain the number of 
muons from a given production interval dz that crosses an arbitrary surface d^A which subtends a solid angle d^fl: 

"^"^^ = d^"^ " = d^ P 
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where ipf^A is the angle between the normal to the surface and the muon direction. On the other hand the projection 
of cPA onto the shower transverse plane is rd(dr — cPAcosipji, where ipA is the angle between the normal to the 
surface and the shower direction. (See Fig. ^Jin Appendix A.) Using this relation we can relate Eq. 01to the number 
of particles per unit area in the transverse plane: 

1 d^NA ^ d^N I cosj^^A ^ d^N 1 ^^^^^ ^ '^"^ « ^ (r.) (6) 

r dzdrdC, dzd'^fl P coatpA dzd^VL P dzdC,dsma sina P 

Where Da{^) denotes the geometrical factor involving these two angles. We stress that Da not only depends on the 
surface orientation with respect to shower axis but also on the direction of the incoming muons. Any detector can be 
regarded as a collection of such surfaces and as a result several such factors Dyi- will have to be considered depending 
on the impact point to obtain the total effective collection area for a given arrival direction. 

We can divide Eq. Elby its integral in z (NArc) to normalize the function to 1. When this is done, using plausible 
functions for the distributions as described in Appendix A, and the relations between r, Z, z and a are used, a number 
of factors cancel out and we obtain the z-distribution of muons arriving at detector A (normalized to 1) which can 
be related to a simple transform of the z distribution: 

1 1 d^NA _ h{z)l^-^ cos a Da{^) 



NatC ^ drdzdC J^' dz h{z) l''^-"' cos a DA{fl)' 

The proposed method relies on the above expression and the geometrical relation between z and t which is described 
in Ref. In that work it was shown that much of the time structure of the muons is due to geometrical effects 

which imply that to each z there corresponds a given arrival time t. As a result we can relate the t and z-distributions: 

1 d^NA _ 1 d^NA dz 
r drd(^dt r drdC,dz dt 

We now define the normalized function g(t) describing the shape of the time distribution through: 

, N 1 d^NA dz 

^A.c5(t) = r:^r:^ ^. (9) 



r drdC^dz 



We can compare the ^-distributions of the muons arriving at ground level given by Eq. |51 to those obtained in 
simulations and agreement is found as will be shown in Section IV. 

The time distribution of the muons is related to the depth distribution of muon production. If we combine Eq. 
and Eq. I^lwe obtain the following relation between h{z) and the time distribution: 



dt 



dz 



1 1 d^NA _ h{z)f~'> cos a DAin) 



9{t)-r =T—z^nrb- TOO (10) 



N^^^rdrdzdC, dzh{z)l^"i cos a Da{^) 



This expression takes into account the fact that from different r we effectively sample the h(z) distribution with an 
extra z-dependence introduced as a overall factor through I and the angles a and C- 

III. RECONSTRUCTION OF THE DEPTH DISTRIBUTION 

In a typical air shower array a number of particle detectors sample the shower front at the Earth's surface. We 
now consider a set of M detectors labeled by a suffix i (from 1 up to M) each with a surface Ai (which becomes Si 
when projected onto the shower plane) and located at position (r, C); in transverse plane coordinates. We calculate 
the time distribution of arriving muons to a detector i by integrating Eq.|Slover the transverse surface ^5*, or, for all 
practical purposes, simply multiplying by the corresponding area Si. Using Eq. |51the time distribution at detector i 
becomes: 



dNA 



dt 



Si 



d^NA, 
dt r dr dC 



{r,0dS^S,NArC9{t). (11) 



The number of muons falling in the detector i can be considered as a finite sample of the continuous arrival time 
distribution probability . Let us assume that we can fill a time histogram with A^^ entries corresponding to the 
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Ni muons detected by detector i. The entries of this histogram can be transformed into a z histogram, using the 
correspondence t ^ z given by [T^ : 



which can, in most cases, be approximated by: 



^(^-^)+A, (12) 



1 r2 



This mapping transforms each time entry (from ~ SiNArC 9{t) ) into a z entry (of Si- f ) and finaUy into 

an entry of the z-distribution of the muons arriving at ground, 

Note that the delay t is the time difference between the arrival time of a given particle and the arrival time of a 
reference plane perpendicular to the shower axis and traveling at the speed of light c, the time-reference plane. We 
have chosen the 0-time origin corresponding to the arrival of the first particle at ground at r = (shower core). If 
the core hits ground at a universal time c^q, the relation between t' and t involves A: 

ct = ct' - ct'o + A. (14) 

Different detectors will give entries to a different time distribution, but they will be converted into samples of a 
unique ^ distribution. As a result the converted entries of available detectors can be combined into a larger sample. 
These entries are naturally weighted by the number of muons detected at each detector. 

In Ref. it was shown that there is an additional source of delay for muons because of their sub-luminal velocities 
t^. The total time delay is obtained adding it to the delay given by Ea.lT^ 

t~tg+t,. (15) 

This delay is energy dependent and it is only dominant over the geometric delay for muons close to shower axis. In 
inclined showers the final muon energy is much larger than the muon mass, mc^ <^ Ei — pal, and therefore: 
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1 (mc^) 



2 cpa 



1 1 



Ei — pal Ei 



(16) 



Since the muon energy is not measured in typical air shower measurements, we cannot account for these effects 
accurately. A solution is simply obtained by eliminating the measurements close to the axis to ensure that the 
kinematical delay can be neglected. This has an impact on statistics. In Fig. ^ we have plotted the factor e{r,z), 
which can be taken as the relative value of the average kinematic delay with respect to the average geometric delay 
(see Appendix B), for ( = 90°, A = and different z. We can see from this graph that the geometric delay dominates 
for distances above 600 m from shower core. 

We can include the kinematic delays on average. We obtain a simple parameterization for the average kinematic 
delay as a function of z and r, (details are given in Appendix B): 

e{r,z)^p,{z)(-y\ (17) 
Vm/ 

If we now subtract the average kinematical delay from the measured delay, instead of Eq. 1121 we obtain: 

1 



2 V ct - c < > 



(ci-c<is>) +A. (18) 



Since it is not possible to obtain z{t) analytically from the previous expression we can use a simple numerical 
approach. We can for instance take zero kinematical delay as a first approximation to obtain z, we then get the average 
kinematical delay and substitute in Eq. 1181 Since the dependence of the p coefBcients on z is mild (logarithmic) the 
procedure converges quickly and one iteration is sufRcient. 
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FIG. 1: The factor e{r,z) which relates to the ratio of kinematic to geometric time delays (see text) versus r for different z: 
1.6, 3.2, 6.3, 13, 25, 50, 100, 200 km, from bottom to top. 



IV. TEST 



One approach to test the method is to simulate showers using a Monte Carlo generator that reproduces the time 
distribution of the signal in a collection of detectors and to apply the method to reconstruct the production distribution 
of the muons h{z). Unfortunately major modifications have to be made in conventional simulation programs that are 
not designed to give the function h{z) to compare with the reconstructed value. We have checked our reconstruction 
method by applying it to showers simulated with the Aires Monte Carlo package 1^. It is straightforward to compare 
the z-distributions of the muons arriving at different locations in the ground. This has been done at different positions 
and agreement is found. This reflects the fact that the time distributions of the muons are well described by the model 
for muon time delays • We have extended the test to combine detectors at different positions and to compare the 
results to the total distribution of the surviving muons dN/dz, which is straightforward to obtain in Aires and in 
most shower generators. In our model this distribution would be obtained by integrating Eq. Qover r and C to cover 
all the ground: 



dN 

dz 



1 d^N 
r drdC^dz 



^drdC, — 



h{z)l^ cos a DA{a) 
Jp°° /i(z)/i~T cosa DA{a) dz 



rdrdC,. 



(19) 



Using simulations we have studied how the total dN/dz distribution at ground relates to the local distributions at 
different r, after integrating in C and r. We first note that to a first approximation the ^-integral is proportional to 
the integrand with A = 0. This is not surprising since A changes sign when integrating over C,. We have found that 
there is an effective value of r (r* ) for local distribution that gives a very good approximation to the overall dN / dz 
distribution. This value is slightly zenith angle dependent and ranges from about 400 m for showers at 0°, up to 
1000 m for horizontal showers at 80° or 1800 m at 86° . This is reasonable because the bulk of the muons arrive to 
ground in a relatively constrained region: for instance, at 0° this region is between 60 m and ~ 1000 m. We can 

and a for a^, = arcsin p- and also consider that to compare with Aires 

i.e. Da = 1. We 



then substitute in Eq. E| I for Z^, — \frl + z^ 
that gives directly the muon position it is not necessary to include a geometric correction factor, 



6 



finally obtain the following approximation 

^ a: h{z)ll-'' cos a^. (20) 

Since in a practical air shower array the detectors are going to be arranged on an unknown and arbitrary pattern 
around the shower axis it is convenient to correct the z-distribution obtained at each detector to a common observation 
distance r^, which approximately reproduces the overall dN/dz as follows: 



dN dt U ^cosa* 1 

— — cx oft) — X — ^ , , ■ (21) 

dz ' dz /i-Tcosa DA{a) 

In general, we must divide by Da{oi) to remove the dependence on the detector geometry if necessary. We note that 

the correction factors /i-t cos a ^^PPi'oach 1 when z increases (i.e. in horizontal showers). Taking this into account 
one can use a single r*- = 400 m for all zenith angles and still obtain relatively good approximations. 

To test the method we have used sets of 500 showers at different zenith angles. A given particle array will be limited 
to a sample of this distribution which is determined by the relative positions of the available detectors. We first take 
the muon output from the simulations and arrange the muons in a time histogram as can be done in an actual shower 
array (we use 25 ns bins). We then apply Eq. ^|to all the simulated muons to calculate z, where we have included 
the geometric corrections of Eq. 1211 and the kinematic corrections as explained in the previous section. Finally an 
r cut is applied. In Figs. |21 and 13 we illustrate the result for a 0° and a 70° showers. The shaded histogram is the 
distribution of all the muon production altitudes compared to that obtained from the reconstruction procedure using 
all the muons which reach the ground with r > Tc- This cut is necessary for the geometric inversion procedure to hold 
accurately. The result indicates that provided the muon time, the shower direction and impact point coordinates are 
known, the reconstruction procedure works well. Figs. [21 and |31 also show the same histogram without the r-cut which 
clearly fails to reproduce the z-distribution obtained in the simulation. This is mostly because of the time accuracy 
of the detectors assumed to be ~ 12.5 ns. 




FIG. 2: Histograms of production distribution for 500 showers of 0° zenith angle and 10^^ eV energy: Light fill: original 
distribution. Unfilled Thiclc Line: Final reconstruction, after all corrections described in the text. Unfilled Thin Line 
Reconstruction with no r cut. 



We have verified that the reconstructed histogram is not sensitive to small changes in the Tf, but clearly these cuts 
in r can have large impact on the statistics. In table we compare the averages of the ratio of z obtained with this 
method to that from the simulation F =< >■ Clearly the method works best for moderately inclined showers 



7 




log^o (z/m) 

FIG. 3: Histograms of production distribution for 500 showers of 70° zenith angle and 10^^ eV energy: Light fill: original 
distribution. Unfilled Thick Line: Final reconstruction, after all corrections described in the text. Unfilled Thin Line 
Reconstruction with no r cut. 



between 60° and 80°. At very low zenith angles, there is an overestimation of the production distance, which could 
be due to an slight overestimation of the energy loss. On the other hand, at very high zenith angles the magnetic field 
effects begin to be important and the time geometric relation underestimates the production distance. Nevertheless, 
in both cases, the precision obtained is quite good. 
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neutrino-like injected at 500 g/cm"^ vertical depth. 

80 I I yes || 1400 | 1.02 

ncutrino-likc injected at 750 g/cm^ vertical depth. 

80 I I yes || 900 | 1.14 



TABLE I: Table of deviation of the reconstruction respect to the real distribution. F =< ^^"^ >. The azimuth angle 4> is 
measured counterclockwise respect to the local magnetic north. For high zenith angles the results obtained with and without 
magnetic field are compared. 

In Fig. ^ we have compared the results of the reconstruction procedure applied to protons and deeply injected 
protons arriving with 80° zenith angle to illustrate how the method can be used to identify deeply interacting inclined 
showers at high zenith, which are natural neutrino candidates. A systematic study of the reconstruction procedure 
and the ability to identify neutrinos under realistic experimental conditions is left for future work. 
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FIG. 4: Histograms of production distribution for 500 showers of 80° zenith angle and 10 eV energy for normal protons and 
for protons injected at 500 r/cm^ of vertical depth to simulate a neutrino interaction (marked as ly-like). The geomagnetic field 
is included. Light fill: original distribution. Unfilled Thick Line: Final reconstruction, after all corrections described in the 
text. 



V. LIMITATIONS OF THE METHOD 

The proposed method gives the reconstruction of the function ^ i.e., the production altitude distribution for 
the surviving muons but clearly the method is only valid with restrictions. Since the method addresses the muon 
time distributions, it is essential to identify the muon contribution in the shower front. Clearly if the muon signal is 
separated at detector level, there are no limitations due to this issue but this is in general not possible for most of the 
detectors used in air shower arrays which typically have complex signals containing a mixture of electron, positron 
and photon signals (the electromagnetic component) and muons. 

Cerenkov detectors, such as the water tanks used in Haverah Park (TBI and in the Auger observatory have some 
advantages in this respect. Since shower muons are penetrating particles the signal they give in Cerenkov detectors is 
basically proportional to their track through the detector, while the electromagnetic component typically gives a signal 
which is simply proportional to the energy carried by photons, electrons and positrons, because it is all absorbed in it. 
As a result the volume of the detector determines the ratio between the muon and electromagnetic signals to a good 
extent. Large detectors give high muon contributions in spite of the muons being a small fraction of all the particles 
in the shower front. Moreover under some circumstances these large sharp pulses could in principle be isolated and 
there are efforts in this line to separate individual muon pulses from the time structure of the Auger tank signals . 

In any case, since the muon lateral distribution is flatter than the electromagnetic contribution p^ . the muons 
eventually dominate for sufficiently large distances to shower axis. As the zenith angle rises the muons dominate 
closer and closer to shower axis. In close to horizontal showers the muons dominate practically always 

Depending on the detector performance there are a number of limitations to the precision which are addressed in 
this subsection. 

It is firstly straightforward to see that the total number of entries of the dN/ dz histogram is the total number of 
muons detected in all the detectors, i.e — '^^^^QNi. If we assume that the ^ distribution has a RMS width cr, 
this means that the position of the mean < z > (which is related to X^ax) can be obtained with no more precision 
than —5= which is an intrinsic statistical limitation. 

A second limitation arises because of the intrinsic time resolution of the detectors used, St, which will limit the 
precision of the muon arrival time. This will translate directly into an uncertainty in the production distance z, Sz, 
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through the map t ^ z. According to Eq. E| 

^.^^^^fl-^U-^. (22) 

z ct \ z J t ^ ' 

This equation can be rewritten to relate — to St substituting ct for the expression given by Eq. 1131 

^ = 2(i^>. (23) 
z z r'^ 

The time resolution of the detector affects the reconstruction precision depending on distance to shower core. As 
we look at the arrival time of muons closer to the shower axis the time delays become smaller and the relative error 
on the z distribution reconstruction diverges. But again this problem can be solved by imposing the cut r > Tc- To 
have — less than a given value e^, and provided that we can find an approximated upper limit for the production 
distance of the muons, z < z^, we obtain the following condition for the cut in r: 



2z„C(5t 



r > r,(C) = . (24) 

1 + J^tan6'cosC 



Here Tc depends on Q because of the asymmetry induced by the term A. If we neglect this term, we obtain a simple 
expression that does not depend on the angle Q: 



2zucSt 

Notice that an r cut can also avoid the regions near the shower axis where the kinematical delay dominates over the 
geometrical, and also the region where the muonic component signal is shadowed by the electromagnetic component. 
In necessary case, the most stringent of the restrictions must be applied. 

For example let us consider an air shower array with time resolution St = 12.5 ns (corresponding to half the sampling 
rate of the Auger detector), located at 1400 m altitude and detecting showers with a zenithal angle of 60°. We can 
easily identify an upper limit for production distance (for instance using simulation), for 60° it is z„ = 31.6 km. 
According to Eq. |55lwe would require that r > 1500 m in order that the resolution on the ^-reconstruction, was 
less that 10% (e^ = 0.1). Fig.[3illustrates the effect showing the contour lines of the precision as a function of r and z 
as given by Eq. 123 The value of Tc must increase as the zenith angle rises because the muons are produced at higher 
z. For high zenith the effect of the neglected ^-dependent term makes the z-reconstruction somewhat worse than the 
approximate expression of Eq. 1251 but enough for our purposes (in necessary case the full expression (Eq. I24() could 
also be used). For a shower with 6 = 86° if we use Cz = 0.1 with the former expression to obtain Tc we actually get a 
resolution — which rises up to 0.17 for the worst case corresponding to ^ = 180°. 



VI. CORRELATION WITH ANGULAR AND CORE POSITION UNCERTAINTIES 



The method has a third intrinsic limitation because to convert the arrival time histogram into a histogram of 
production distances using Eq. lAllI the incoming direction and the position of shower axis must be known, so that 
the appropriate values of r can be introduced. In a realistic case these will only be known to a given precision and 
further uncertainties will arise because of the correlations between both core position and direction uncertainties 
with the distance distribution obtained. The study of these effects with simulations shows that the reconstructed 
distributions have a minimum in width when the true shower direction and impact position are used to reconstruct 
the depth distribution. This adds a interest to the method since it can be used in principle as a further check of the 
reconstructed directions and impact points. 

To study these correlations we explore the stability of the reconstruction to shifts in the core positions and angular 
directions. Unfortunately the computing time necessary to test such stability can be very large if simulations are used 
in the same way they were used to test the method in the previous section. We will use instead the results of the muon 
time delay model of Ref. to get distributions of the arrival time for the shower muons from simulated showers. In 
an attempt to be closer to experimental conditions we assume an array of particle detectors and calculate the number 
of muons that crosses each of them. We choose 10 x 1.2 m (area x height) detectors in a hexagonal grid, separated 
1500 m corresponding to the Auger surface detector. This limits the statistics of the reconstructed distribution, dN/dz, 
in a realistic way. Fig. shows an example of the statistics that could be obtained for a 10^^ eV shower with 9 = 60° 
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FIG. 6: Distribution for a reconstructed lO'^^ eV proton shower of ^ = 60° and (f> = 90° (histograms) compared with the original 
distribution when we have finite sampling produced by a finite number of detectors. 



and — 90°. ( The azimuth angle was measured counterclockwise respect to East direction.) We first recalculate 
the depth distributions assuming that the incidence direction has been misreconstructed by A0) with respect 
to the actual arrival direction chosen for the simulation. This procedure was repeated for different shifts in angular 
space within an interval of 4° x 4° . For each angular shift both the mean and RMS width of the z-distributions in 
log basis were calculated. In Figs. |3 and |H1 representative results showing the mean and RMS width for an example 
of a 10^^ eV proton shower of 60° zenith are shown as a function of Ad (A0) for fixed (p (6), Also the RMS width is 
shown in a two dimensional plot Ac/)) in Fig. |H1 It is worth remarking that the mean value of z is quite stable 
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FIG. 7: Effects of reconstruction direction shifts for a 10^^ eV proton shower oi 6 — 60° and (j> = 90° with an r cut r > 1500 m. 
Mean of \ogiQ{z/m) as a function of A9 fixing A<j) = (dashed line) and as a function of Aij) fixing A9 = (continuous hue). 
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FIG. 8: Effects of reconstruction direction shifts in the width of the z-distribution for a 10^^ eV proton shower oi = 60° and 
(j) = 90° with an r cut r > 1500 m. Left RMS width of logj^Q(2;/m) as a function of angular shifts {A9,A<j)). Right RMS of 
logig{z / m) as a function of fixing A(j) = (dashed line) and as a function of A(j) fixing A6 = (continuous line). 



to shifts in azimuthal direction A0, whereas there is a very slight rise of the reconstructed z when increasing the 
zenith angle 9. The behavior on this angle in not symmetrical since A introduces an asymmetry between early and 
late regions. On the other hand the RMS width of the distributions seems to have a local minimum when the correct 
arrival direction is used. In a typical air shower array the arrival direction is obtained using the relative arrival times 
of the signals at different locations. The observed minimum of the RMS width of the z-distribution suggests that 
this method could be also used to either reconstruct shower direction independently or, more likely, to check that 
the reconstruction obtained through conventional methods is consistent with the arrival time of the muons at large 
distances from shower core, on an individual shower basis. 

A completely analogous method was followed to study the core position and z-reconstruction correlations. The 
reconstructed impact points were shifted by (Ax, Ay) in the ground plane with respect to the core of the simulated 
shower over a grid covering a rectangle of 1000 x 1000 m. The means and RMS widths for shifts in both x and y 
positions are shown in Fig. |51 The plots display important discontinuities. They are of statistical nature because the 
total number of muons in the detector is small and as the core position position is shifted individual detectors are 
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FIG. 9; Effects of impact point shifts on the mean (left panel) and RMS (right panel) of the ^-distribution as a function of Ax 
fixing Ay — (dashed line) and as a function of Ay fixing Aa; = (continuous line) for a 10^^ eV proton shower at 6 — 60° 
and — 90° with an r cut r > 1500 m. 



rejected or accepted because of the r cut. The detectors close to the cut are those that have most muons and including 
or not including them affects the z-distribution. These discontinuities are also present in some circumstances for an- 
gular shifts because the relative position of the detectors also change but clearly the changes of angular reconstruction 
modify the distances to shower axis at a second order level. These discontinuities can become smoother by increasing 
statistics, for instance relaxing the r-cut. 

Notice that the mean value of z is again quite stable to shifts in core position. This is not difficult to understand 
since each core location corresponds to a new time-reference plane which is only slightly shifted along the shower axis 
with respect to the planes obtained for other core positions. The effect is due to the curvature of the front and is 
therefore a second order effect. The RMS width of the distributions also displays a local minimum when the correct 
impact point is used. This is also not surprising given that approximately z oc . Differences in the z reconstruction 
arise through the modification of the relative position of the tank with respect to shower core. It can be easily seen 
that when a tank gets a closer position (r decreases) as a result of the shift, the tanks placed in the opposite C will 
get to a further one (r increases). This suggests that the width of the reconstructed distribution should have a local 
minimum when the correct impact point is considered. An example is given relaxing the r-cut to 500 m in Figs. lTUI 

In a typical air shower array the core position is obtained by the relative amplitude of the signals at different locations 
either through a fit or by some other means. As for shifts in angular direction the minimum of the RMS width of the 
z-distribution suggests that this method could be used either to reconstruct shower impact points independently or 
to check that the core position reconstruction obtained through conventional methods is consistent with the arrival 
time of the muons at large distances from shower core, on an individual shower basis. 



VII. SUMMARY 



We have developed a method that has the potential of reconstructing the production altitude for the muons in 
inclined cosmic ray showers based on the time distribution of the muon signals in the detectors of an extensive air 
shower array. This method requires knowledge of the arrival times of muons in the detectors of the air shower array 
and it can be applied provided that a cut is made in distance to shower axis, r > rc- Since the muon signal dominates 
at high r it can be also used when the detectors cannot separate the muon signal provided that the r cut is chosen so 
that the muon signal dominates. 

The method relies mainly on geometrical arguments and there are minor effects introduced through the kinematical 
delay of the muons which have little effect at large distances from shower axis. The model does not rely on any as- 
sumption about the interaction model for hadrons. Different models would give rise to different kinematic corrections, 
but the effect is small. Although we have assumed proton showers to explore the viability of this method the method 
can be also used for heavier nuclei, and similar results would be obtained in that case. The necessary cut introduces 
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FIG. 10: Effects of impact point shifts on the RMS width of the ^-distribution for a 10^^ eV proton at 6 — 60° and (j> = 90° 
with an relaxed r cut r > 500 m . Left RMS of logj^Q(2:/m) as a function of a core shift (Aa;, Ay). Right RMS of \ogj^Q{z/m) 
as a function of Aa; (dashed line) fixing Ay = 0, and as a function of Ay fixing Ax = (continuous line). 



limitations because of statistics. We have checked that our method correctly reproduces the depth distribution of 
muon production using sets of simulated showers with AIRES to a degree of accuracy that is zenith angle dependent. 
The method works best in the 60° — 80° region and it is fairly stable with respect to misreconstruction of the shower 
core and the incoming direction, in what concerns the mean of the distribution. The RMS width of the production 
distance distribution however is sensitive to both the reconstructed impact point and arrival direction. The width 
displays a minimum when the correct impact point and arrival direction are used in the reconstruction procedure. 

This work represents a new approach to studying extensive air showers. It will add information concerning the 
individual development of air showers and can be used to check the reconstructed impact point and arrival directions. 
The reconstruction of depth development in inclined showers can also have important implications in improving the 
potential of air shower arrays to detect neutrinos. 
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APPENDIX A: MODELLING THE DISTRIBUTION OF SURVIVING MUONS 



At ground level both muon energies and muon number are reduced because of energy loss and decay. As a first 
approximation, assuming a constant energy loss per unit depth dE/dx = —a along an uniform atmospheric density 
p, both these effects can be easily accounted for. After traveling a distance I a flux of muons 4>q of energy reduces 
by an energy dependent factor to: 



Ei — pal 
E, 



(Al) 



The last factor takes into account both energy loss and decay in flight. The muon mean lifetime, r, enters through 
the exponent k = mc^/{pacT) ~ 0.8. We can correct the energy of the produced muons given by Eq. ^with such 
factor to obtain the distribution in Ei, pt, z and ^ after traveling a distance I: 



d'^N(l) 1 



Ei — pal 



E, 



(A2) 



Here I is the distance traveled by the muon, which enters in Eq. lAll given by P = {z — A)^ + r^, where r is the 
distance to shower axis at the end of the muon trajectory, and z — A is the distance travelled by the muon measured 
along the shower axis. The correction A — r tan 6* cos C relates the distances measured along the axis between the 
intercepts of shower axis and the muon trajectory with ground level and depends on the zenith angle, 0, as well as on 
the muon impact point coordinates in the transverse plane, r and C. 

For a muon to reach ground level there is a minimal production energy given by Ei > mc^ + pal. Typical values of 
pal greatly exceed mc^, particularly for inclined showers. After traveling a distance I the transverse distance is simply 
r = Zsina. Performing the change the coordinates from pt to sin a in Eg. IA2l we obtain: 



d^N{l) 



1 



dzdQdEidsma 2n 



E,. 



1 - 



pal 

'e~ 



E^ 

c 



(A3) 



Correlations between the ground variables appear naturally because of energy loss and decay. 

We can now introduce the following parameterizations for fi{Ei) and f2{pt) which were shown in to give good 
approximations to the muon time distributions: 



fliE^) 



7 



1 



E, 



where 7 ~ 2.6 and Q ~ 0.17 GeV/c. Eg. IA3l now becomes 



d-^N 



dz d( dEi dsin a 



^foh-l) 

2-K 



hiz) 



E, 



-7+1 



Ei sin a 
,2n2 



C2Q2 



Ei sin a 



pal 

'e~ 



We now integrate the distribution in Ei for fixed z and a to obtain: 



d^N 



dz d(^ dsin a 



d'^N 



c^+nai d(dEi dsin a 



■dE, 



2tt 



Qc 



in which /(/, sin a) is a dimensionless integral in the variable x = ^'^^q " : 

I{l,sma) ~ / x^'^^'^ 1 eyip{—x)dx, 



(A4) 
(A5) 

(A6) 

(A7) 

(A8) 



15 



with Xo{l, ct) — ^ sin a and y{a) = sin a. We note that y xq when / mc^/pa ~ 500 m and in that case the 
integral can be approximated replacing its lower limit by xq. Since z, the distance of muon production, is relatively 
large, particularly for inclined showers, this approximation is adequate for most circumstances. Then it is easily seen 
that the integral depends only on the combination I sin a = r, i.e. on the transverse distance. We can then replace 
I {I, sin a) by /(r). 

In terms of the differential solid angle (Pil — —d(^d cos a the number of muons arriving at ground level coming from 
production distance z becomes: 




FIG. 11: Illustration of the relation between the ground surface, the transverse plane (shower plane) and the shower axis (big 
arrow) . The small arrows represent a detector surface direction normal to the detector surface, the dashed lines are parallel to 
shower axis and the continuous lines are the muon trajectories. The expressions for number of muons (INa which go through 
the surface element cPA (labeled 1) are described in the text. 

Using the relations between z, I and a, we can substitute Eg . IA9I into Eq. Inland integrate in z to obtain: 
^ _ld^NA AAo(7-l) fmc^y-' I{r) f ^ cos a DAjn) sin^-' a ^ 

NatC is the number of muons going through d^A (whose projection onto the transverse plane is rdrd() and has a 
non-trivial geometrical dependence through the Da{^) factor. This factor can be both greater or less than 1 and 
becomes 1 as the angle a tends to zero, that is in the limit of small r. In that case NatC simply becomes the muon 
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surface density in the transverse plane. Da{^) can be regarded as a geometrical correction which accounts for the fact 
that the muons do not travel parallel to the shower axis and the fact that the detector planes are not perpendicular 
to the shower axis not to the muon fluxes (See Fig. 111(1 . 

It can be shown that when z and r are fixed Da depends on the angle ^ because the angles a and 4'fj.A differ. As a 
result the factor Da is responsible for part of the asymmetry in the muon signal. For instance a horizontal surface will 
have larger efficiency for collecting early arriving muons because ipp^A is smaller than for late arriving muons. Note 
that Da{^) would be exactly 1 for an "ideal" detector of spherical shape. In that case these differences disappear. 

The expression obtained in Ea. llUI relates the time distribution to a transform of the depth distribution h{z) which 
effectively accounts for muon decay in flight through I. If we are interested in obtaining the production distribution 
h(z) formally we can rewrite Ea. llUI as: 

(if r°° 

h{z) ^ g(t)—l-^+'^ cos-^ a D^\n) hiz)!^"^ cos a DA{n)dz. (All) 
dz Jq 

This expression in principle allows us to obtain h{z) from the time distribution of the arriving muons at a given point 
on the ground, g{t). On its own it is not very useful because typically a single detector in an air shower array does 
not collect sufficient statistics to sample h{z) reliably, in a practical situation one must combine the results of several 
detectors. Since h(z) is normalized to 1, the unknown factor which is given by the integral on the RHS of the equation 
acts as an effective weight to be given to each detector. In a first attempt it is possible to ignore it. For inclined 
showers the weights to be applied for the relevant detectors are expected to be quite similar and the approximation 
works well. More sophisticated approaches could be deviced for instance using Eg. IAll| with a trial h{z) function in 
an iterative process to sample h{z). However since h(z) is not directly available from the simulation program used we 
will not need to calculate h(z) and we have instead compared our results to the z-distribution of the surviving muons. 



APPENDIX B: PARAMETERIZATION OF KINEMATICAL DELAYS 

The mean kinematical time delay can be obtained by applying the method developed in Ref. p^ . summarized in 
EqElto the distributions discussed in this article: 

ft d'*Af ^jT- 

, J ^ dzdC,dEidr ^ ^-^^^\ 

J dzdCdEidr^'- 

After some manipulation, using the results of the models in Appendix A a simple expression can be obtained for it: 

lr^(m^\^ Iylx„ [1 -^Y~^ exp {-x) dx i ^2 



2c I \cQ J /^_^^^ a;-'r+2 [l - f] exp {-x) dx 2c I 

with xo{l, ct) = ^ sina and y{a) = sina. 

In the last equality of the above expressions we have introduced the dimensionless factor e(r, z — A). We note that 
for z — A ^ r, for instance in inclined showers, it gives the ratio of the average kinematical delay to the geometrical 
delay at a given position ~ e(r, z — A). This indicates the regions where the geometric delay dominates, which 

depend on production distance. For practical purposes we have parameterized e(r, z) as follows: 

e{r,z)=poiz)(-Y\ (B3) 
Vm/ 

with 

logiopo(^) = -0.6085+ 1.955 logio(z/m) - 0.3299 logJo(z/m) + 0.0186 logJo(z/m), (B4) 
logioPi = -1.176. (B5) 



